Appunti di Statistica Spaziale

1 carica librerie

library(sf)

## Warning: il pacchetto 'sf' è stato creato con R versione 4.4.2

## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

require(ggplot2)

## Caricamento del pacchetto richiesto: ggplot2

library(tidygeocoder)

## Warning: il pacchetto 'tidygeocoder' è stato creato con R versione 4.4.2

library(dplyr)  

## Warning: il pacchetto 'dplyr' è stato creato con R versione 4.4.2

## 
## Caricamento pacchetto: 'dplyr'

## I seguenti oggetti sono mascherati da 'package:stats':
## 
##     filter, lag

## I seguenti oggetti sono mascherati da 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spatstat)

## Warning: il pacchetto 'spatstat' è stato creato con R versione 4.4.3

## Caricamento del pacchetto richiesto: spatstat.data

## Warning: il pacchetto 'spatstat.data' è stato creato con R versione 4.4.3

## Caricamento del pacchetto richiesto: spatstat.univar

## Warning: il pacchetto 'spatstat.univar' è stato creato con R versione 4.4.3

## spatstat.univar 3.1-2

## Caricamento del pacchetto richiesto: spatstat.geom

## Warning: il pacchetto 'spatstat.geom' è stato creato con R versione 4.4.3

## spatstat.geom 3.3-5

## Caricamento del pacchetto richiesto: spatstat.random

## Warning: il pacchetto 'spatstat.random' è stato creato con R versione 4.4.3

## spatstat.random 3.3-2

## Caricamento del pacchetto richiesto: spatstat.explore

## Warning: il pacchetto 'spatstat.explore' è stato creato con R versione 4.4.3

## Caricamento del pacchetto richiesto: nlme

## 
## Caricamento pacchetto: 'nlme'

## Il seguente oggetto è mascherato da 'package:dplyr':
## 
##     collapse

## spatstat.explore 3.3-4

## Caricamento del pacchetto richiesto: spatstat.model

## Warning: il pacchetto 'spatstat.model' è stato creato con R versione 4.4.3

## Caricamento del pacchetto richiesto: rpart

## spatstat.model 3.3-4

## Caricamento del pacchetto richiesto: spatstat.linnet

## Warning: il pacchetto 'spatstat.linnet' è stato creato con R versione 4.4.3

## spatstat.linnet 3.2-5

## 
## spatstat 3.3-1 
## For an introduction to spatstat, type 'beginner'

library(scatterplot3d)
library(rgl)

## Warning: il pacchetto 'rgl' è stato creato con R versione 4.4.2

library(stars)

## Warning: il pacchetto 'stars' è stato creato con R versione 4.4.3

## Caricamento del pacchetto richiesto: abind

library(gstat)

## Warning: il pacchetto 'gstat' è stato creato con R versione 4.4.3

## 
## Caricamento pacchetto: 'gstat'

## Il seguente oggetto è mascherato da 'package:spatstat.explore':
## 
##     idw

library(plotly)

## 
## Caricamento pacchetto: 'plotly'

## Il seguente oggetto è mascherato da 'package:ggplot2':
## 
##     last_plot

## Il seguente oggetto è mascherato da 'package:stats':
## 
##     filter

## Il seguente oggetto è mascherato da 'package:graphics':
## 
##     layout

library(sp)

2 carica i path delle cartelle

base = 'C:/Users/feder/OneDrive/Documenti/Fede/Unimib/statistica spaziale'
paths = paste0(base, '/lab',1:5)

3 funzioni utili

3.0.1 sf

st_read : carica un file .shp come oggetto sf

st_write : esporta un oggetto sf in un file .shp

st_set_crs : setta il sistema di riferimento (crs=coordinate reference system) di un oggetto sf (che è anche un dataframe con una geometry) I codici più comuni sono 3003 (Gauss-Boaga fuso Ovest) e 4326 (WGS84).

st_transform: permette di cambiare il sistema di riferimento di un oggetto sf attenzione: con st_set_crs si importa l’oggetto con un crs che si deve sapere, mentre con st_transform si cambia il crs di un oggetto già importato

st_crs: permette di vedere il sistema di riferimento di un oggetto sf ($input)

st_geometry: permette di estrarre la geometria di un oggetto sf st_boundary : calcola il perimetro del rettangolo in cui è circoscritto un oggetto sf geom_sf : permette di visualizzare un oggetto sf in ggplot

st_centroid : dato un oggetto sf (che è una lista di poligoni), calcola il centroide di ogni poligono. Se l’oggetto sf è un poligono, il centroide è il punto medio del poligono.

st_intersection : dati due oggetti sf, restituisce l’oggetto sf dove il secondo è contenuto nel primo

st_join : unisce due oggetti sf (o sp) in base a una colonna comune

st_as_sf: prende un dataframe, e lo trasforma in un oggetto sf, chiedendo le coordinate x e y (coords = c(“x”, “y”)) e il sistema di riferimento (crs = 3003 per gauss-boaga fuso ovest, 32632 per wgs84 zona 32)

st_bbox: si ottengono i 4 estremi laterali di un oggetto sf (le coordinate xmin ymin xmax ymax delle altezze del rettangolo in cui è circoscritto l’oggetto poly)

3.0.2 tidygeocoder

geo: prende un dataframe con indirizzi e restituisce le coordinate geografiche i valori richiesti sono: - street - city - postalcode - country Mentre lat e lon sono due stringhe arbitrarie (i nomi da dare alle colonne della latitudine e longitudine del dataframe di output) method = ‘osm’ di default, e indica l’API service

3.0.3 sp

merge: unisce due oggetti sf (o sp) in base a una colonna comune

3.0.4 spatstat

as.ppp: crea un oggetto ppp (point pattern) a partire da un dataframe a due colonne e da una finestra W che è un oggetto sf

ppp: crea un oggetto ppp (point pattern) specificando le coordinate x e y, la finestra window (as.owin(oggetto sf)) e il marker (se numerico il plot farà una scala di dimensioni, se factor il plot farà una legenda di diversi tipi di punti. Puoi modificare i colori, le dimensioni e le forme dei punti con i parametri col, cex e pch)

ppp(df$Est, df$Nord, window=as.owin(lomb.poly), marks=df$Ettari)
ppp(df$Est, df$Nord, window=as.owin(lomb.poly), marks=as.factor(df$Ettari > median(df$Ettari)) 

disc: crea una finestra circolare per il processo. Richiede radius = scalare e center = c(x,y)

3.0.5 esempio di sf

setwd(paths[1])


inc0 = incendi = read.table("Incendi_2003.csv", header=T,sep=";")
head(inc0)

##   Ettari    Nord     Est
## 1    9.0 5076900 1604750
## 2    0.1 5098880 1518920
## 3    2.0 5090650 1508800
## 4    0.1 5083300 1491200
## 5    0.2 5083300 1491300
## 6    2.0 5075000 1563520

inc1 = st_as_sf(inc0, coords = c("Est", "Nord"), crs = 3003)

lomb.poly<- st_read("lombardia.shp", quiet=TRUE) 
lomb.poly= st_set_crs(lomb.poly, 3003) #non cambia la classe, solo le coordinate
class(lomb.poly)

## [1] "sf"         "data.frame"

ggplot(data = st_boundary(lomb.poly)) +  #plotta il oundary estratto a oggetto st
  geom_sf()

ggplot(data = inc1) +
  geom_sf()

ggplot() +
  geom_sf(data = st_boundary(lomb.poly))+ 
  geom_sf(data = inc1,  size = 2.5) + 
  theme_bw()  +
  ggspatial::annotation_north_arrow(which_north = "true") +
  ggspatial::annotation_scale(location="br")+
  ggtitle('Incendi in Lombardia 2003')

##senza ggplot

#in gauss boaga fuso ovest
plot(st_boundary(st_geometry(lomb.poly)),lwd=4,col='violetred', 
     main='Lombardia Gauss Boaga senza ggplot')
points(incendi$Est,incendi$Nord)

#in wgs84 zona 32
lomb.polyUTM<-st_read("Lombardia_UTMWGS84.shp",quiet=TRUE) 
lomb.polyUTM= st_set_crs(lomb.polyUTM, 32632)  #UTM-WGS84 zona 32

plot(st_boundary(st_geometry(lomb.polyUTM)),lwd=4,col='red', 
     main='Lombardia UTM WGS84 senza ggplot e senza punti')
points(incendi$Est,incendi$Nord, col='green') #i punti non ci sono!! perché sono in un crs diverso

#per vedere i punti in un altro crs, bisogna trasformarli
plot(st_boundary(st_geometry(lomb.polyUTM)),lwd=4,col='blue',
     main='Lombardia UTM WGS84 senza ggplot e con i punti')
inc2 = st_transform(inc1, crs = 32632) #trasforma il crs dell'oggetto
points(inc2$geometry, col='green') #ora i punti ci sono

print('st_bbox')

## [1] "st_bbox"

st_bbox(lomb.polyUTM)

##      xmin      ymin      xmax      ymax 
##  460659.1 4947434.7  691622.0 5165401.8

plot(st_boundary(st_geometry(lomb.polyUTM)),lwd=4,col=2,
     xlim = c(460659,1691467 ),   ##allarga l'asse x
     ylim = st_bbox(lomb.polyUTM)[c(2, 4)],  #lascia uguale l'asse y
     main='ecco perché prima non vedevi i punti'
)
points(incendi$Est,incendi$Nord)

3.0.6 esempio di sf 2

##estraggo ora il poligono della Lombardia dal poligono dell'Italia
##lo faccio con filter, perché un oggetto sf è anche un dataframe

setwd(paths[1])

reg.poly<-st_read("Reg01012023_g_WGS84.shp",quiet=TRUE);
print('st_crs(reg.poly)$input')

## [1] "st_crs(reg.poly)$input"

st_crs(reg.poly)$input 

## [1] "WGS 84 / UTM zone 32N"

plot(st_boundary(st_geometry(reg.poly)), main='italia')

head(as.data.frame(reg.poly)) #vedi i nomi delle colonne)

##   COD_RIP COD_REG               DEN_REG Shape_Leng  Shape_Area
## 1       1       1              Piemonte  1236799.8 25393881930
## 2       1       2         Valle d'Aosta   310968.1  3258837561
## 3       1       3             Lombardia  1410222.9 23862315752
## 4       2       4   Trentino-Alto Adige   800893.7 13607548167
## 5       2       5                Veneto  1054587.0 18343552568
## 6       2       6 Friuli Venezia Giulia   670820.7  7934115959
##                         geometry
## 1 MULTIPOLYGON (((457749.5 51...
## 2 MULTIPOLYGON (((390652.6 50...
## 3 MULTIPOLYGON (((485536.4 49...
## 4 MULTIPOLYGON (((743267.7 52...
## 5 MULTIPOLYGON (((768124 5175...
## 6 MULTIPOLYGON (((872344.5 50...

lomb = reg.poly %>% filter(COD_REG == 3) #la lombardia nel poligono ha codice 3
#ma se facevo  lomb=reg.poly %>% filter(DEN_REG == 'Lombardia') era uguale
print('st_crs(lomb)$input')

## [1] "st_crs(lomb)$input"

st_crs(lomb)$input

## [1] "WGS 84 / UTM zone 32N"

plot(st_boundary(st_geometry(lomb)), main='lombardia')

# esportare il nuovo shape
  #st_write(lomb, "lomb_mia.shp") #salvo il mio shape file
  #check : si vede che riesco ad aprire e leggere lo shape file
  #a<-st_read("lomb_mia.shp",quiet=TRUE) 
  #a= st_set_crs(a, 32632)
  #ggplot(data = st_boundary(a)) + 
  #  geom_sf(); rm(a)

3.0.7 esempio di sf 3

  centr = st_centroid(reg.poly)  # estrazione centroidi

## Warning: st_centroid assumes attributes are constant over geometries

#il centroide dev'essere dentro il poligono. 
# Per poligoni diagonali, come la Liguria, 
# il centroide calcolato fuori come media grezza
# dal poligono esce dal poligono stesso, per cui va aggiustato idoneamente
  
  
  
  ggplot(data = st_boundary(reg.poly)) + 
    geom_sf(data = centr,  size = 2.5)+  #qui si fa il grosso
    ggspatial::annotation_north_arrow(which_north = "true") +
    ggspatial::annotation_scale(location="br") +   
    geom_sf_text(data = centr,    
                 aes(label = DEN_REG),     
                 size = 2.5,     
                 color = "blue" ,     
                 nudge_x = -2000,     
                 nudge_y = 25000) +  
    geom_sf()  #si vede che il centroide della Liguria è in mare:(

  # osservare: i centroidi ereditano le info su CRS:  st_crs(centr)
  st_crs(centr)$input

## [1] "WGS 84 / UTM zone 32N"

  st_crs(reg.poly)$input

## [1] "WGS 84 / UTM zone 32N"

3.0.8 esempio di tidygeocoder

setwd(paths[1])

data_firm=dd=read.csv("datafirm.csv", sep=";", header = T)
data_firm = data_firm[1:5,]
position<-geo(street = data_firm$address, 
              postalcode = data_firm$cap, 
              city=data_firm$city,
              country = data_firm$country, 
              method = 'osm', lat = myla , long = mylo)

## Passing 5 addresses to the Nominatim single address geocoder

## Query completed in: 5.1 seconds

poscoords = st_as_sf(position, coords = c("mylo","myla"), crs=4326) 
#poscoords = st_transform(poscoords, crs = 3003) #trasforma il crs dell'oggetto
st_crs(poscoords)$input

## [1] "EPSG:4326"

reg.poly<-st_read("Reg01012023_g_WGS84.shp",quiet=TRUE)
st_crs(reg.poly)$input #controlla il crs dell'oggetto: è diverso da quello sopra

## [1] "WGS 84 / UTM zone 32N"

reg.poly= st_transform(reg.poly, 4326)
st_crs(reg.poly)$input #ora è corretto

## [1] "EPSG:4326"

plot(st_boundary(st_geometry(reg.poly)),col='violetred', 
     main='Posizioni delle aziende senza ggplot')
points(poscoords$geometry, col='green', lwd=2)

ggplot() +
  geom_sf(data = st_boundary(st_geometry(reg.poly)), col='violetred') +
  geom_sf(data = poscoords, col='green', lwd=2) + 
  theme_bw() +
  ggspatial::annotation_north_arrow(which_north = "true") +
  ggspatial::annotation_scale(location="br")+
  ggtitle('Posizioni delle aziende con ggplot')

## Scale on map varies by more than 10%, scale bar may be inaccurate

3.0.9 esempio sf 4

setwd(paths[1])

OMI<-st_read("MI_Zone_OMI corretto_region.shp",quiet=TRUE); st_crs(OMI)

## Coordinate Reference System:
##   User input: Monte Mario / Italy zone 1 
##   wkt:
## PROJCRS["Monte Mario / Italy zone 1",
##     BASEGEOGCRS["Monte Mario",
##         DATUM["Monte Mario",
##             ELLIPSOID["International 1924",6378388,297,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4265]],
##     CONVERSION["Italy zone 1",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",9,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.9996,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",1500000,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Italy - onshore and offshore - west of 12°E."],
##         BBOX[36.53,5.93,47.04,12]],
##     ID["EPSG",3003]]

  OMI= st_set_crs(OMI, 3003)  # Gauss-Boaga
    #plot(st_boundary(st_geometry(OMI)))
  MM<-st_read("MM_FERMATE.shp", quiet=TRUE);st_crs(MM)

## Coordinate Reference System: NA

  MM= st_set_crs(MM,3003)  # Gausso-Boaga
  
  plot(st_boundary(st_geometry(OMI)))

  ggplot() + geom_sf(data=MM)

  ggplot() + geom_sf(data = st_geometry(OMI)) + geom_sf(data=MM)

  #ggplot(data = st_boundary(OMI)) + geom_sf() + geom_sf(data=MM)
  
  MMinMI <- st_intersection(OMI, MM)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

  ggplot(data = st_boundary(OMI)) + geom_sf() + geom_sf(data=MMinMI) 

  U7<-geo(street ="Via Bicocca degli Arcimboldi", postalcode = , city="Milano",
          country = "Italy", method = 'osm')

## Passing 1 address to the Nominatim single address geocoder

## Query completed in: 1 seconds

  U7 <- sf::st_as_sf(U7, coords = c("long","lat"),crs=4326)
    st_crs(U7)$input

## [1] "EPSG:4326"

  #utm_proj <- st_crs(3003)  
  U7 <- st_transform(U7, crs = 3003) 

  size = 1000  #espresso in metri
  buffer    = st_buffer(U7, dist = size) 
  buffer_sf = st_as_sf(buffer)  
  
  
  ggplot() + 
    geom_sf(data = st_boundary(OMI)) + 
    geom_sf(data=MMinMI) + 
    geom_sf(data=U7, size=2, col='blue') +
    geom_sf(data=st_boundary(buffer_sf), size=2, col='red', shape=3) 

    #con st_boundary ho l'interno vuoto
    
  
# Esercizio: Ripetere l'operazione di buffering attorno all'U7 4 volte per un raggio dell'intorno 
#            pari a 0.75, 1, 1.25, 1.5 chilometri. Riportare i buffer con colori diversi sulla stessa mappa
  
  sizes = c(0.75, 1, 1.25, 1.5) * 1000 #espresso in metri
  
  mappa = ggplot() + 
    geom_sf(data = st_boundary(OMI)) + 
    geom_sf(data=MMinMI) + 
    geom_sf(data=U7, size=2, col='blue')
  
  colori = rainbow(length(sizes))
  for (i in 1:length(sizes)){
    size = sizes[i]
    buffer    = st_buffer(U7, dist = size)
    buffer_sf = st_as_sf(buffer)
    mappa = mappa + geom_sf(data=st_boundary(buffer_sf), size=2, col=colori[i], shape=3)
  }
  print(mappa)

#come contare il numero di fermate in ogni zona OMI  
a <- st_join(OMI, MM, left = T)
stops <- a %>% count(ZONA) #da dplyr, come table ma su una colonna del df


ggplot(data=stops) + 
    geom_sf(aes(fill=n, geometry=geometry))+
    scale_fill_continuous(low="yellow", high="red", na.value="gray")+
    labs(title = "metro stops in Milan")+ 
    ggspatial::annotation_north_arrow(which_north = "true") +
    ggspatial::annotation_scale(location="br") +
    geom_sf(data=MM, color = "blue", size = 1.5)

## nuovo esempio  
  
aus.poly = st_read("austrianuts3.shp", quiet=TRUE)
aus.poly = st_set_crs(aus.poly, 4326)  # WGS84
head(aus.poly)  

## Simple feature collection with 6 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 14.45184 ymin: 46.83054 xmax: 17.15986 ymax: 48.38746
## Geodetic CRS:  WGS 84
##      ID                    NAME                       geometry
## 1 AT111        Mittelburgenland POLYGON ((16.28007 47.45503...
## 2 AT112          Nordburgenland POLYGON ((16.37145 47.64218...
## 3 AT113           Südburgenland POLYGON ((15.99618 46.83518...
## 4 AT121 Mostviertel-Eisenwurzen POLYGON ((14.73689 47.74818...
## 5 AT122    Niederösterreich-Süd POLYGON ((15.21628 47.79579...
## 6 AT123            Sankt Pölten POLYGON ((15.33461 47.89, 1...

distkm <- read.table("dist.txt", header=T)
head(distkm)

##      id id_num  dist pop.dens
## 1 AT111    111 18806   100.16
## 2 AT112    112 18573   102.25
## 3 AT113    113 21311   117.83
## 4 AT121    121 18507   128.16
## 5 AT122    122 17819   223.18
## 6 AT123    123 16791   173.40

dati.sp = sp::merge(aus.poly, distkm, by.y = "id", by.x ="ID", all.x=T)

# draw a map
 ggplot(data=dati.sp)+
  geom_sf(aes(fill=dist, geometry=geometry))+
   scale_fill_continuous(low="yellow", high="red", na.value="gray")+
    labs(title = "Car use - average distances driven by car")+ 
    ggspatial::annotation_north_arrow(which_north = "true") +
    ggspatial::annotation_scale(location="br") 

3.0.10 esempio di spatstat

  setwd(paths[2])

  lomb.poly<-st_read("lombardia.shp",quiet=TRUE)  #leggi il file shape con gli allegati
  lomb.poly= st_set_crs(lomb.poly, 3003)  #3003 codice gauss boaga

  incendi <- read.table("Incendi_2003.csv", header=T,sep=";");  
  str(incendi)      

## 'data.frame':    380 obs. of  3 variables:
##  $ Ettari: num  9 0.1 2 0.1 0.2 2 1.15 0.05 0.02 3 ...
##  $ Nord  : int  5076900 5098880 5090650 5083300 5083300 5075000 5075200 5114790 5116640 4952750 ...
##  $ Est   : int  1604750 1518920 1508800 1491200 1491300 1563520 1638900 1601120 1603410 1524300 ...

  ppp<-incendi[,c("Est","Nord")]        #escludo il marker ettari
  #crea un oggetto ppp (point pattern) a partire da un oggetto sf
  ppp0 = as.ppp(ppp, W=lomb.poly);ppp0 

## Warning: data contain duplicated points

## Planar point pattern: 380 points
## window: polygonal boundary
## enclosing rectangle: [1460708, 1691467] x [4947383, 5165460] units

  print('info sull oggetto ppp')

## [1] "info sull oggetto ppp"

  class(ppp0)

## [1] "ppp"

  summary(ppp0)

## Planar point pattern:  380 points
## Average intensity 1.592947e-08 points per square unit
## 
## *Pattern contains duplicated points*
## 
## Coordinates are integers
## i.e. rounded to the nearest unit
## 
## Window: polygonal boundary
## 4 separate polygons (no holes)
##            vertices        area relative.area
## polygon 1     12245 23852200000      1.00e+00
## polygon 2        99     2655710      1.11e-04
## polygon 3        31      177156      7.43e-06
## polygon 4        27      171204      7.18e-06
## enclosing rectangle: [1460708, 1691467] x [4947383, 5165460] units
##                      (230800 x 218100 units)
## Window area = 23855200000 square units
## Fraction of frame area: 0.474

  # plot PP
  plot(ppp0,cex=0.5,main="Incendi in Lombardia nel 2003")

# PPS marcati con dimensione incendio
  df<-incendi[,c("Est","Nord","Ettari")]
  ppp0=ppp(df$Est, df$Nord, window=as.owin(lomb.poly), marks=df$Ettari);ppp0

## Warning: data contain duplicated points

## Marked planar point pattern: 380 points
## marks are numeric, of storage type  'double'
## window: polygonal boundary
## enclosing rectangle: [1460708, 1691467] x [4947383, 5165460] units

  #as.owin prende un oggetto geometrico e lo forza in modo tale da poterlo mostrare come un 
  #oggetto raffigurabile nella finestra del point pattern
  plot(ppp0,
       main="Incendi in Lombardia nel 2003 differenziati per estensione")

mark = df$Ettari > median(df$Ettari)
  
# come factor binary: cambia l'estetica
ppp1 = ppp(df$Est,df$Nord, window=as.owin(lomb.poly),
           marks = as.factor(ifelse(mark, 'hell', 'candle')))

## Warning: data contain duplicated points

ppp1

## Marked planar point pattern: 380 points
## Multitype, with levels = candle, hell 
## window: polygonal boundary
## enclosing rectangle: [1460708, 1691467] x [4947383, 5165460] units

plot(ppp1,
     main="Incendi in Lombardia nel 2003 differenziati per estensione",
     cols = c('darkred', 'orange')[as.numeric(mark)+1],
     pch = c(17, 4)[as.numeric(mark)+1],
     lwd=2
)

3.0.11 spatstat: test sull’ipotesi di CSR di un ppp

setwd(paths[2])

  d<-read.table("WLR300.txt",header=T)  # lettura dataset esterno in formato testo
  str(d)

## 'data.frame':    66 obs. of  2 variables:
##  $ X: int  -6064 -2785 -5827 -5978 -11416 -3581 -1636 27189 14599 30251 ...
##  $ Y: int  1663 7358 6091 7326 -10350 15263 17727 965 24172 -5110 ...

  W <- disc(96950, c(0,0))          # definisce una finestra circolare per il processo
  #disc è una funzione di spatstat che crea oggetti di classe finestra
  pp<-ppp(d$X,d$Y,window=W); summary(pp)

## Planar point pattern:  66 points
## Average intensity 2.236005e-09 points per square unit
## 
## Coordinates are integers
## i.e. rounded to the nearest unit
## 
## Window: polygonal boundary
## single connected closed polygon with 128 vertices
## enclosing rectangle: [-96950, 96950] x [-96950, 96950] units
##                      (193900 x 193900 units)
## Window area = 29516900000 square units
## Fraction of frame area: 0.785

  plot(pp,cex=0.6,main="locazioni eventi",cex.main =1.2,lwd=3)

  qx<-quadratcount(pp, 4, 4)    ### test basato sul quadrat counts:  tabella 4x4
    plot(qx) #si partiziona la regione, e si verifica se in ogni cella cade una certa osservazione